Transformation method for diffusion spectrum imaging using large deformation diffeomorphic metric mapping

ABSTRACT

A transformation method for diffusion spectrum imaging includes: receiving an original DSI dataset and a template DSI dataset; computing an energy function; computing, for each time point, first-order and second-order derivatives of the energy function with respect to velocity fields in an image space and in a q-space; computing, for each time point, the velocity fields in the image space and in the q-space based upon the first-order and second-order derivatives; performing integration on the velocity fields over time to obtain a deformation field; and generating a transformed DSI dataset according to the deformation field.

CROSS-REFERENCE TO RELATED APPLICATION

This application claims priority to Taiwanese Application No. 101129000,filed on Aug. 10, 2012.

BACKGROUND OF THE INVENTION

1. Field of the Invention

The invention relates to an image processing method, and moreparticularly to an image processing method for magnetic resonanceimaging.

2. Description of the Related Art

Diffusion MRI is a non-invasive imaging method suitable for evaluatingfiber orientation of a specific region and revealing the underlyingwhite matter structure of the human brain. One of the approaches is thediffusion tensor imaging (DTI), where the water diffusion is modelled asGaussian distribution. Practically, a DTI dataset consists of sixdiffusion weighted (DW) images by applying the diffusion-sensitivegradients in six non-colinear directions, and one null image with noapplication of diffusion-sensitive gradient. The water diffusion isencoded in the signal intensity of the DW images. Therefore, withappropriate processing on the signal intensity of the DTI dataset, thediffusion tensor in each voxel can be estimated. The diffusion tensorcould be represented by a symmetric 3-by-3 matrix, where the principaleigenvector of the matrix is usually assumed to coincide with theunderlying fiber orientation. The fiber orientation map can be furtherprocessed to reconstruct the fiber pathways.

However, the Gaussian assumption limits DTI to detect at most one fiberorientation in each voxel. Consequently, in regions with crossingfibers, it is difficult for DTI to resolve the fiber orientations andwould lead to inaccurate estimation of the anisotropy index.

The crossing fiber problem could be resolved through estimating thediffusion orientation distribution function (ODF). The diffusion ODFcould be estimated by the high angular resolution diffusion image(HARDI) methods, such as q-ball imaging (QBI), or by a grid samplingscheme, which is also called diffusion spectrum imaging (DSI). All ofthe methods do not impose any diffusion models.

In neuroimage studies, it is usually required to transform the brainimages to a common space, the so-called template space, to perform theanalyses (e.g., statistical comparisons between healthy and patientgroups). It is known in the art to use a linear or non-linear method totransform a three-dimensional (3D) brain image to the template space.However, the conventional 3D transformation methods can only deal withscalar images. For an appropriate transformation on diffusion imagessuch as DTI, QBI or DSI, not only the anatomical structures need to beregistered, but the diffusion profiles require to be aligned.

SUMMARY OF THE INVENTION

Therefore, an object of the present invention is to provide atransformation method that may completely transform image spatialinformation and diffusion information for diffusion spectrum imaging(DSI) datasets.

According to the present invention, a transformation method for DSIcomprises:

a) receiving an original DSI dataset and a template DSI dataset;

b) using a processor to compute an energy function E;

c) using the processor to obtain, for each time point, a first-orderderivative and a second-order derivative of the energy function Ecomputed in step b) with respect to velocity field in an image space anda first-order derivative and a second-order derivative of the energyfunction E computed in step b) with respect to velocity field in aq-space;

d) using the processor to compute, for each time point, the velocityfield in the image space and the velocity field in the q-space basedupon the first-order and second-order derivatives obtained in step c);

e) using the processor to perform integration on the velocity fieldsobtained in step d) over time to obtain a deformation field, which mapsthe original DSI dataset to the template DSI dataset; and

f) using the processor to generate a transformed DSI dataset accordingto the deformation field obtained in step e).

BRIEF DESCRIPTION OF THE DRAWINGS

Other features and advantages of the present invention will becomeapparent in the following detailed description of the preferredembodiment with reference to the accompanying drawings, of which:

FIG. 1 is a flow chart illustrating steps of a preferred embodiment ofthe transformation method for diffusion spectrum imaging (DSI) accordingto the present invention;

FIG. 2 is an image of an original DSI of an example;

FIG. 3 is an image of the template DSI of the example; and

FIG. 4 is an image transformed using the preferred embodiment of thetransformation method.

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENT

Referring to FIG. 1, the preferred embodiment of the transformationmethod for diffusion spectrum imaging (DSI) according to this inventionis implemented by a processor of a computer after the computer loads aproprietary program stored in a storage medium. The transformationmethod comprises the following steps:

Step S11: Receiving an original DSI dataset (such as the image shown inFIG. 2) and a template DSI dataset (such as the image shown in FIG. 3).A DSI image set often includes hundreds of images, and FIGS. 2 and 3 areobtained by summarizing the hundreds of the images to a generalizedfractional anisotropy (GFA) map.

Step 12: Using a processor of the computer to compute an energy functionE, which is:

$\begin{matrix}{E = {{E_{1} + E_{2}} = {\underset{\underset{E_{1}}{︸}}{\frac{1}{2}{\int_{0}^{1}{{v_{t}}^{2}\ {\mathbb{d}t}}}} + \underset{\underset{E_{2}}{︸}}{\frac{1}{2\sigma^{2}}{\int{\left( {{\left( {WI}_{0} \right) \cdot g_{10}} - {WI}_{1}} \right)^{2}{\mathbb{d}q}{\mathbb{d}x}}}}}}} & (1)\end{matrix}$where E₁ is an energy of the transformation path, E₂ is an energyrepresenting data-matching, x is a three-dimensional coordinate in animage space, and q is a three-dimensional coordinate in a q-space. v_(t)is velocity field at time t, σ is a parameter controlling the relativecontribution between E₁ and E₂, W is a predetermined weighting function,I₀ is the template DSI dataset, and I₁ is the original DSI dataset. Theinterpretation of g_(ab) is that for a particle located at (x,q) at timet=a, g_(ab)(x,q) is a position of the particle at time t=b. Thus, g₁₀could be considered as a deformation mapping function of WI₀ when t=1.The weighting function W(q) is to compensate low signal intensity inregions with high |q| values.

Step 13: Using the processor to obtain, for each time point, afirst-order derivative and a second-order derivative of the energyfunction E computed in step 12 with respect to velocity field in theimage space and a first-order derivative and a second-order derivativeof the energy function E computed in step 12 with respect to velocityfield in the q-space.

A large deformation diffeomorphic metric mapping (LDDMM), which isproposed by Beg, M. F. et al. in 2005, is used in this preferredembodiment. The transformation process is assumed to behave like liquidflow, so that there is an associated velocity field at each time point.In other words, the velocity field is a function of time.

In this embodiment, the first-order derivative of the energy function Ewith respect to the velocity field in the image space is computed usingthe following equation (2):

$\begin{matrix}{\frac{\partial E}{\partial v_{x,t}} = {v_{x,t} - {K\left\lbrack {\frac{1}{\sigma^{2}}{\int{{{Dg}_{t,1}}{\nabla_{x}{J_{t}^{0}\left( {J_{t}^{0} - J_{t}^{1}} \right)}}{\mathbb{d}q}}}} \right\rbrack}}} & (2)\end{matrix}$

The first-order derivative of the energy function E with respect to thevelocity field in the q-space is computed using the following equation(3):

$\begin{matrix}{\frac{\partial E}{\partial v_{q,t}} = {v_{q,t} - {K\left\lbrack {\frac{1}{\sigma^{2}}{{Dg}_{t,1}}{\nabla_{q}{J_{t}^{0}\left( {J_{t}^{0} - J_{t}^{1}} \right)}}} \right\rbrack}}} & (3)\end{matrix}$where v_(x,t) is the velocity field in the image space at time t,v_(q,t) is the velocity field in the q-space at time t, K is a smoothingoperator characterizing the smoothness of the velocity field, J_(t)⁰=(WI₀)∘g_(t0), and J_(t) ¹=(WI₁)∘g_(t1).

The second-order derivative of the energy function E with respect to thevelocity field in the image space is computed using the followingequation (4):

$\begin{matrix}{\frac{\partial^{2}E}{\partial v_{x,t}^{2}} = {K\left\lbrack {\frac{1}{\sigma^{2}}{\int{{{Dg}_{t,1}}\left( {\nabla_{x}J_{t}^{0}} \right)\left( {\nabla_{x}J_{t}^{0}} \right)^{T}{\mathbb{d}q}}}} \right\rbrack}} & (4)\end{matrix}$

The second-order derivative of the energy function E with respect to thevelocity field in the q-space is computed using the following equation(5):

$\begin{matrix}{\frac{\partial^{2}E}{\partial v_{q,t}^{2}} = {K\left\lbrack {\frac{1}{\sigma^{2}}{{Dg}_{t,1}}\left( {\nabla_{q}J_{t}^{0}} \right)\left( {\nabla_{q}J_{t}^{0}} \right)^{T}} \right\rbrack}} & (5)\end{matrix}$

Step 14: Using the processor to compute, for each time point, thevelocity field in the image space and the velocity field in the q-spacebased upon the first-order and second-order derivatives obtained in step13.

In this embodiment, the first-order and second-order derivativesobtained in step 13 are applied to a Levenberg-Marquardt (LM) algorithmto accelerate calculation as in the following equations (6) and (7):

$\begin{matrix}{v_{x,t}^{({n + 1})} = {v_{x,t}^{(n)} - {K\left\lbrack {\int{\left( {{\frac{1}{\sigma^{2}}{{Dg}_{t,1}}\left( {\nabla_{x}J_{t}^{0}} \right)\left( {\nabla_{x}J_{t}^{0}} \right)^{T}} + {\delta\; I}} \right)^{- 1}\left( {{K^{- 1}v_{x,t}^{(n)}} - {\frac{1}{\sigma^{2}}{{Dg}_{t,1}}{\nabla_{x}{J_{t}^{0}\left( {J_{t}^{0} - J_{t}^{1}} \right)}}}} \right){\mathbb{d}q}}} \right\rbrack}}} & (6) \\{v_{q,t}^{({n + 1})} = {v_{q,t}^{(n)} - {K\left\lbrack {\left( {{\frac{1}{\sigma^{2}}{{Dg}_{t,1}}\left( {\nabla_{q}J_{t}^{0}} \right)\left( {\nabla_{q}J_{t}^{0}} \right)^{T}} + {\delta\; I}} \right)^{- 1}\left( {{K^{- 1}v_{q,t}^{(n)}} - {\frac{1}{\sigma^{2}}{{Dg}_{t,1}}{\nabla_{q}{J_{t}^{0}\left( {J_{t}^{0} - J_{t}^{1}} \right)}}}} \right)} \right\rbrack}}} & (7)\end{matrix}$where K⁻¹ is an inverse operator of K. Calculation may converge afterabout 5 iterations (n=5), and the velocity field in the image space andthe velocity field in the q-space at time t are thus obtained.

Step 15: Using the processor to perform integration on the velocityfields obtained in step 14 over time to obtain a deformation field whichmaps the original DSI dataset to the template DSI dataset. Thedeformation field is obtained according to the following equation (8):

$\begin{matrix}{\frac{\mathbb{d}g_{0t}}{\mathbb{d}t} = {v_{t} \circ g_{0t}}} & (8)\end{matrix}$

where ∘ denotes function composition. Here the deformation fieldg_(0t)=(g_(x,0t), g_(q,0t)), where g_(x,0t), which is an image spacecomponent of g_(0t), is a function of x, and g_(q,0t), which is aq-space component of g_(0t), is a function of both x and q. Similarly,v_(t)=(v_(x,t), v_(q,t)), where v_(x,t), which is an image spacecomponent of v_(t), is a function of x, and v_(q,t), which is a q-spacecomponent of v_(t), is a function of both x and q.

Step 16: Using the processor to generate a transformed DSI datasetaccording to the deformation field obtained in step 15. The transformedDSI dataset is shown as the image in FIG. 4.

It should be noted that the processors used in each step of this methodmay be the same or different.

To sum up, the transformation method of this invention uses not only thethree-dimensional information in the image space, but also thethree-dimensional diffusion information of the DSI dataset, in a totalof six dimensions. Therefore, when the original DSI dataset isregistered to the template space, details like diffusion information inthe image space and in the q-space may be completely transformed, so asto enhance precision of the subsequent comparison and application.

While the present invention has been described in connection with whatis considered the most practical and preferred embodiment, it isunderstood that this invention is not limited to the disclosedembodiment but is intended to cover various arrangements included withinthe spirit and scope of the broadest interpretation so as to encompassall such modifications and equivalent arrangements.

What is claimed is:
 1. A transformation method for diffusion spectrumimaging (DSI) comprising: a) receiving an original DSI dataset and atemplate DSI dataset; b) using a processor to compute an energy functionE; c) using the processor to obtain, for each time point, a first-orderderivative and a second-order derivative of the energy function Ecomputed in step b) with respect to velocity field in an image space anda first-order derivative and a second-order derivative of the energyfunction E computed in step b) with respect to velocity field in aq-space; d) using the processor to compute, for each time point, thevelocity field in the image space and the velocity field in the q-spacebased upon the first-order and second-order derivatives obtained in stepc); e) using the processor to perform integration on the velocity fieldsobtained in step d) overtime to obtain a deformation field, which mapsthe original DSI dataset to the template DSI dataset, wherein thedeformation field is calculated according to:${\frac{\mathbb{d}g_{0t}}{\mathbb{d}t} = {v_{t} \circ g_{0t}}};$ andwherein the deformation field g_(0t)=(g_(x,0t), g_(q,0t)), where x is athree-dimensional coordinate in the image space, q is athree-dimensional coordinate in the q-space, g_(x,ot), which is an imagespace component of g_(0t), is a function of x, and g_(q,0t), which is aq-space component of g_(0t), is a function of both x and q; and thevelocity field v_(t)=(v_(x,t), v_(q,t)) where v_(x,t), which is an imagespace component of v_(t), is a function of x, and v_(q,t), which is aq-space component of v_(t), is a function of both x and q; and f) usingthe processor to generate a transformed DSI dataset according to thedeformation field obtained in step e).
 2. The transformation method asclaimed in claim 1, wherein, in step d), the velocity field for eachtime point is calculated using an iterative Levenberg-Marquardtalgorithm.
 3. The transformation method as claimed in claim 1, wherein,in step b), the energy function is: $\begin{matrix}{E = {{E_{1} + E_{2}} = {\underset{\underset{E_{1}}{︸}}{\frac{1}{2}{\int_{0}^{1}{{v_{t}}^{2}\ {\mathbb{d}t}}}} + \underset{\underset{E_{2}}{︸}}{\frac{1}{2\sigma^{2}}{\int{\left( {{\left( {WI}_{0} \right) \circ g_{10}} - {WI}_{1}} \right)^{2}{\mathbb{d}q}{\mathbb{d}x}}}}}}} & (1)\end{matrix}$ where E₁ is an energy of the transformation path, E₂ is anenergy representing data-matching, v_(t) is the velocity field at timet, σ is a parameter controlling the relative contribution between E₁ andE₂, W is a predetermined weighting function, I₀ is the template DSIdataset, I₁ is the original DSI dataset, g₁₀ is a deformation mappingfunction of WI₀ when t=1, x is a three-dimensional coordinate in theimage space, and q is a three-dimensional coordinate in the q-space. 4.A transformation method for diffusion spectrum imaging (DSI) comprising:a) receiving an original DSI dataset and a template DSI dataset; b)using a processor to compute an energy function E: $\begin{matrix}{E = {{E_{1} + E_{2}} = {\underset{\underset{E_{1}}{︸}}{\frac{1}{2}{\int_{0}^{1}{{v_{t}}^{2}\ {\mathbb{d}t}}}} + \underset{\underset{E_{2}}{︸}}{\frac{1}{2\sigma^{2}}{\int{\left( {{\left( {WI}_{0} \right) \circ g_{10}} - {WI}_{1}} \right)^{2}{\mathbb{d}q}{\mathbb{d}x}}}}}}} & (1)\end{matrix}$ where E₁ is an energy of the transformation path, E₂ is anenergy representing data-matching, v_(t) is the velocity field at timet, σ is a parameter controlling the relative contribution between E₁ andE₂, W is a predetermined weighting function, I₀ is the template DSIdataset, I₁ is the original DSI dataset, g₁₀ is a deformation mappingfunction of WI₀ when t=1, x is a three-dimensional coordinate in theimage space, and q is a three-dimensional coordinate in the q-space; c)using the processor to obtain, for each time point, a first-orderderivative and a second-order derivative of the energy function Ecomputed in step b) with respect to velocity field in an image space anda first-order derivative and a second-order derivative of the energyfunction E computed in step b) with respect to velocity field in aq-space; d) using the processor to compute, for each time point, thevelocity field in the image space and the velocity field in the q-spacebased upon the first-order and second-order derivatives obtained in stepc); e) using the processor to perform integration on the velocity fieldsobtained in step d) overtime to obtain a deformation field, which mapsthe original DSI dataset to the template DSI dataset; f) using theprocessor to generate a transformed DSI dataset according to thedeformation field obtained in step e).
 5. The transformation method asclaimed in claim 4, wherein, in step e), the deformation field iscalculated according to:${\frac{\mathbb{d}g_{0t}}{\mathbb{d}t} = {v_{t} \circ g_{0t}}};$ wherev_(t) is velocity field at time t, g_(0t) is the deformation field attime t, and ° denotes function composition.
 6. The transformation methodas claimed in claim 4, wherein, in step d), the velocity field for eachtime point is calculated using an iterative Levenberg-Marquardtalgorithm.